Importing data

designMatrix <- ReadDataFrameFromTsv("../../design/all_samples_short_names.tsv.csv")
## ../../design/all_samples_short_names.tsv.csv read from disk!
exonmatrix <- ReadDataFrameFromTsv(file.name.path="../../data/NEW_PFC_exonMatrix.txt", 
                            row.names.col=NULL)
## ../../data/NEW_PFC_exonMatrix.txt read from disk!
counts <- exonmatrix[,c(7:dim(exonmatrix)[2])]

exonDesignMatrix <- cbind(colnames(counts), designMatrix)
exonDesignMatrix$samples <- rownames(exonDesignMatrix)
exonDesignMatrix$exonsamples <- colnames(counts)
rownames(exonDesignMatrix) <- exonDesignMatrix$exonsamples

exonDesignMatrix <- exonDesignMatrix[-c(grep("RS2",exonDesignMatrix$condition)),]
exonDesignMatrix<- exonDesignMatrix[-c(grep("HC7",exonDesignMatrix$condition)),]


exonmatrix <- exonmatrix[, c(1:6, which(colnames(exonmatrix) %in% rownames(exonDesignMatrix)))]
colnames(exonmatrix)
##  [1] "GeneID"  "Chr"     "Start"   "End"     "Strand"  "Length"  "S3HC5_1"
##  [8] "S3HC5_2" "S3HC5_3" "S3HC5_4" "S3HC5_6" "S3SD5_1" "S3SD5_2" "S3SD5_3"
## [15] "S3SD5_4" "S3SD5_6" "WTHC5_1" "WTHC5_2" "WTHC5_3" "WTHC5_4" "WTHC5_6"
## [22] "WTSD5_1" "WTSD5_2" "WTSD5_3" "WTSD5_4" "WTSD5_6"
annot <- exonmatrix[,c(1:6)]
counts <- exonmatrix[,c(7:dim(exonmatrix)[2])]

y.all <- DGEList(counts=counts, genes=annot)

y.all <- y.all[, colnames(counts)]
summary(is.na(y.all$genes))
##    GeneID           Chr            Start            End         
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:222996    FALSE:222996    FALSE:222996    FALSE:222996   
##    Strand          Length       
##  Mode :logical   Mode :logical  
##  FALSE:222996    FALSE:222996
dge <- y.all

loading annotation file

## annotation downloaded from ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/

ncbi.annotation <- read.delim("../../data/annotation/Mus_musculus.gene_info", 
                    comment.char = "%", header=TRUE, stringsAsFactors=FALSE)

map <- match(dge$genes$GeneID, ncbi.annotation$GeneID)
# map <- map[!is.na(map)]
# print("ERROR: NA ON CHROMOSOME")
# dge$genes$newChr <- ncbi.annotation$chromosome[map] 
dge$genes$Chr <- gsub(pattern="chr", replacement="", x=dge$genes$Chr)
dge$genes$Symbol <- ncbi.annotation$Symbol[map]
dge$genes$Strand <- NULL
# dge$genes
dge.annot <- dge
summary(is.na(dge.annot$genes))
##    GeneID           Chr            Start            End         
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:222996    FALSE:222996    FALSE:222996    FALSE:222996   
##                                                                 
##    Length          Symbol       
##  Mode :logical   Mode :logical  
##  FALSE:222996    FALSE:220491   
##                  TRUE :2505

filtering data

library(dplyr)
library(tidyr)

long_data <- gather(exonmatrix, key=condition, value=expression, 
                rownames(exonDesignMatrix))

long_data %>%
    group_by(GeneID, condition) %>%
    summarize(total = sum(expression)) %>%
    group_by(GeneID) %>%
    summarize(pass = sum(total > 20) >= 5) %>%
    filter(pass) %>%
    pull(GeneID) -> keep_ids

keep <- which(exonmatrix$GeneID %in% keep_ids)
# length(keep)
# dim(exonmatrix)

dge.annot <- dge.annot[keep, , keep.lib.sizes=FALSE]

Normalizations

TMM

dge.norm <- calcNormFactors(dge.annot)
dge.counts <- estimateCommonDisp(dge.norm)

PlotPCAPlotlyFunction(counts.data.frame=dge.counts$pseudo.counts, 
            design.matrix=exonDesignMatrix, shapeColname="condition", 
            colorColname="genotype", plotly.flag=TRUE, show.plot.flag=TRUE)
## [1] TRUE

loading negative controls

## negative controls
library(readxl)

sd.ctrls <- read_excel(path="../../data/controls/Additional File 4 full list of BMC genomics SD&RS2.xlsx", sheet=1)
sd.ctrls <- sd.ctrls[order(sd.ctrls$adj.P.Val),]

sd.neg.ctrls <- sd.ctrls[sd.ctrls$adj.P.Val > 0.7, ]

sd.neg.ctrls <- sd.neg.ctrls$`MGI Symbol`
sd.neg.ctrls <- sd.neg.ctrls[-which(is.na(sd.neg.ctrls))]

tab <- table(dge.counts$genes$Symbol)
single <- names(which(tab==1))
negcon <- which(dge.counts$genes$Symbol %in% intersect(sd.neg.ctrls, single))
# length(negcon)

exonDesignMatrix$gcondition <- droplevels(exonDesignMatrix$gcondition)

ruving

library(RUVSeq)
groups <- makeGroups(exonDesignMatrix$gcondition)
ruvedSExprData <- RUVs(as.matrix(round(dge.counts$pseudo.counts)), cIdx=negcon,
                       scIdx=groups, k=5)

normExprData <- ruvedSExprData$normalizedCounts

PlotPCAPlotlyFunction(counts.data.frame=normExprData, 
            design.matrix=exonDesignMatrix, shapeColname="condition", 
            colorColname="genotype", plotly.flag=TRUE, show.plot.flag=TRUE,
            prefix.plot="TMM+RUVs")
## [1] TRUE
plotRLE(normExprData, outline=FALSE)

edgering with ruv

design <- model.matrix(~0+exonDesignMatrix$gcondition + ruvedSExprData$W)

colnames(design) <- c(levels(exonDesignMatrix$gcondition), paste0("W", 1:5))

dge.disp.des <- estimateDisp(dge.norm, design, robust=TRUE)
# dge.disp.des$pseudo.counts
# dge.counts$common.dispersion
# plotBCV(dge.disp.des)
fit <- glmFit(dge.disp.des, design)

KOSD5 - WTSD5

Differential expression

contr <- limma::makeContrasts(contrasts="KOSD5 - WTSD5", levels=design)
qlf <- glmLRT(fit, contrast=contr)
topTags(qlf, p.value=0.05, n=20)
## Coefficient:  1*KOSD5 -1*WTSD5 
##           GeneID Chr     Start       End Length  Symbol      logFC
## 72408      58234  15  89547626  89549882   2257  Shank3 -8.7534513
## 72409      58234  15  89557735  89560261   2527  Shank3  0.7104138
## 72404      58234  15  89533374  89533454     81  Shank3  1.0999900
## 72657      12805  15  92220447  92220479     33   Cntn1  1.8655127
## 72397      58234  15  89521103  89521373    271  Shank3  0.6915776
## 73190     634104  15  98220808  98221056    249 Olfr287 -6.4774186
## 72406      58234  15  89543100  89543230    131  Shank3  0.7775398
## 72402      58234  15  89532095  89532219    125  Shank3  0.6528003
## 73189     634104  15  98210004  98210123    120 Olfr287 -4.0643742
## 72389      58234  15  89500215  89500418    204  Shank3  0.7091774
## 72414      11434  15  89573832  89574585    754     Acr  1.1831918
## 72393      58234  15  89503540  89503707    168  Shank3  0.6926872
## 72654      12805  15  92128128  92128335    208   Cntn1 -6.0994072
## 73188     634104  15  98206971  98208377   1407 Olfr287 -3.2129539
## 72394      58234  15  89503797  89503913    117  Shank3  0.7546330
## 49819  100039826  13  12259813  12260273    461    <NA>  0.4544742
## 131933     73284   3 137625966 137628332   2367  Ddit4l  0.3714430
## 72392      58234  15  89503303  89503454    152  Shank3  0.6588258
## 72399      58234  15  89525160  89525273    114  Shank3  0.6481363
## 170987     68498   6 127949793 127953031   3239 Tspan11  0.5560605
##            logCPM         LR       PValue          FDR
## 72408   5.7438646 5912.87004 0.000000e+00 0.000000e+00
## 72409   7.2956758  173.38308 1.349924e-39 1.172004e-34
## 72404   2.7280258  130.33269 3.465448e-30 2.005802e-25
## 72657   0.0711896   70.37475 4.904355e-17 2.128981e-12
## 72397   3.3468886   63.03623 2.029395e-15 7.047684e-11
## 73190  -2.0498492   61.69349 4.013055e-15 1.161378e-10
## 72406   2.4094219   59.17273 1.444212e-14 3.582470e-10
## 72402   3.0971448   56.45919 5.737707e-14 1.245369e-09
## 73189  -1.8661725   52.51822 4.262748e-13 8.224263e-09
## 72389   2.5887533   49.98179 1.551797e-12 2.694541e-08
## 72414   1.1416585   49.43803 2.047383e-12 3.231886e-08
## 72393   2.4987683   47.97280 4.321734e-12 5.971279e-08
## 72654  -2.1793780   47.90642 4.470550e-12 5.971279e-08
## 73188  -0.1910068   41.63178 1.101873e-10 1.366637e-06
## 72394   1.7228002   39.25924 3.711082e-10 4.295948e-06
## 49819   3.6101228   33.71623 6.376655e-09 6.920264e-05
## 131933  4.8442894   33.58932 6.806534e-09 6.952274e-05
## 72392   1.9592242   33.41622 7.440104e-09 7.177220e-05
## 72399   1.7188706   28.52627 9.243592e-08 8.439784e-04
## 170987  2.2821088   28.42877 9.721013e-08 8.439784e-04

Alternative Splicing

sp <- diffSpliceDGE(fit, contrast=contr, geneid="GeneID", exonid="Start")
## Total number of exons:  173640 
## Total number of genes:  15630 
## Number of genes with 1 exon:  796 
## Mean number of exons in a gene:  11 
## Max number of exons in a gene:  314
topSpliceDGE(sp, test="Simes", FDR=0.05)
##       GeneID Chr Symbol NExons      P.Value          FDR
## 72409  58234  15 Shank3     22 0.000000e+00 0.000000e+00
## 70774  27367  15   Rpl3     10 1.763053e-20 1.307657e-16
## 72678  12805  15  Cntn1     26 1.976317e-13 9.772230e-10
## 73447  22143  15 Tuba1b      4 8.382376e-09 3.108604e-05

exon plot

plotSpliceDGE(sp, geneid="Shank3", genecol="Symbol")

plotSpliceDGE(sp, geneid="Cntn1", genecol="Symbol")

plotSpliceDGE(sp, geneid="Rpl3", genecol="Symbol")

plotSpliceDGE(sp, geneid="Tuba1b", genecol="Symbol")

plotSpliceDGE(sp, geneid="Fabp7", genecol="Symbol")

plotSpliceDGE(sp, geneid="Hdac7", genecol="Symbol")

plotSpliceDGE(sp, geneid="Per3", genecol="Symbol")

plotSpliceDGE(sp, geneid="Kalrn", genecol="Symbol")

plotSpliceDGE(sp, geneid="Sgk1", genecol="Symbol")

plotSpliceDGE(sp, geneid="Dido1", genecol="Symbol")

edgering without ruv (QLfit)

design <- model.matrix(~0+exonDesignMatrix$gcondition)

colnames(design) <- levels(exonDesignMatrix$gcondition)

dge.disp.des <- estimateDisp(dge.norm, design, robust=TRUE)

fit <- glmFit(dge.disp.des, design)#glmQLFit(dge.disp.des, design)

KOSD5 - WTSD5

Differential expression (QLtest)

contr <- limma::makeContrasts(contrasts="KOSD5 - WTSD5", levels=design)
qlf <- glmLRT(fit, contrast=contr)#glmQLFTest(fit, contrast=contr)
topTags(qlf, n=20, p.value=0.05)
## Coefficient:  1*KOSD5 -1*WTSD5 
##        GeneID Chr     Start       End Length  Symbol      logFC
## 72408   58234  15  89547626  89549882   2257  Shank3 -8.8022894
## 72409   58234  15  89557735  89560261   2527  Shank3  0.7182385
## 73190  634104  15  98220808  98221056    249 Olfr287 -6.6280910
## 73188  634104  15  98206971  98208377   1407 Olfr287 -3.2926207
## 72654   12805  15  92128128  92128335    208   Cntn1 -6.2038428
## 72404   58234  15  89533374  89533454     81  Shank3  1.0362329
## 73189  634104  15  98210004  98210123    120 Olfr287 -4.2429463
## 72414   11434  15  89573832  89574585    754     Acr  1.2101703
## 72406   58234  15  89543100  89543230    131  Shank3  0.8154886
## 72657   12805  15  92220447  92220479     33   Cntn1  1.7057011
## 72402   58234  15  89532095  89532219    125  Shank3  0.6449115
## 73883  109901  15 100681143 100681288    146   Cela1 -4.0288217
## 72389   58234  15  89500215  89500418    204  Shank3  0.6897990
## 124624  26570   3  50364936  50372366   7431 Slc7a11  0.5835040
##             logCPM         LR       PValue          FDR
## 72408   5.74935846 2581.80484 0.000000e+00 0.000000e+00
## 72409   7.29565025   70.88405 3.788489e-17 3.289166e-12
## 73190  -2.04039759   59.77164 1.065272e-14 5.972940e-10
## 73188  -0.17532367   59.26803 1.375937e-14 5.972940e-10
## 72654  -2.18106438   53.62246 2.429664e-13 8.437737e-09
## 72404   2.73493738   51.50843 7.129051e-13 2.063147e-08
## 73189  -1.85552807   45.81991 1.296394e-11 3.215799e-07
## 72414   1.12615618   36.07217 1.901434e-09 4.127062e-05
## 72406   2.40620387   34.66236 3.921403e-09 7.565693e-05
## 72657   0.08011166   32.94057 9.501940e-09 1.649917e-04
## 72402   3.09506433   28.48043 9.465025e-08 1.494097e-03
## 73883  -2.50180153   24.24342 8.489681e-07 1.228457e-02
## 72389   2.59030798   23.97110 9.779280e-07 1.306211e-02
## 124624  5.80058770   21.32691 3.872572e-06 4.803096e-02

Alternative Splicing

sp1 <- diffSpliceDGE(fit, contrast=contr, geneid="GeneID", exonid="Start")
## Total number of exons:  173640 
## Total number of genes:  15630 
## Number of genes with 1 exon:  796 
## Mean number of exons in a gene:  11 
## Max number of exons in a gene:  314
topSpliceDGE(sp1, test="Simes", FDR=0.05)
##       GeneID Chr Symbol NExons      P.Value          FDR
## 72409  58234  15 Shank3     22 0.000000e+00 0.000000e+00
## 70774  27367  15   Rpl3     10 8.409808e-22 6.237554e-18
## 72678  12805  15  Cntn1     26 1.275235e-12 6.305610e-09
## 73447  22143  15 Tuba1b      4 9.852492e-08 3.653797e-04
WriteDataFrameAsTsv(data.frame.to.save=topSpliceDGE(sp1, test="Simes", FDR=1, 
            number=Inf), file.name.path="KOSD5_WTSD5_diff_Splice_edger")

exon plot

plotSpliceDGE(sp1, geneid="Shank3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Cntn1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Rpl3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Tuba1b", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Fabp7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Hdac7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Per3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Kalrn", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Sgk1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Dido1", genecol="Symbol")

WTSD5 - WTHC5

Differential expression (QLtest)

contr <- limma::makeContrasts(contrasts="WTSD5 - WTHC5", levels=design)
qlf <- glmLRT(fit, contrast=contr)#glmQLFTest(fit, contrast=contr)
topTags(qlf, n=20, p.value=0.05)
## Coefficient:  -1*WTHC5 1*WTSD5 
##           GeneID Chr     Start       End Length  Symbol     logFC
## 89251      57435  17  56103417  56105501   2085   Plin4  4.026792
## 89248      57435  17  56100591  56102330   1740   Plin4  4.033977
## 89252      57435  17  56106096  56107374   1279   Plin4  3.524487
## 219969     14605   X 140542846 140543185    340 Tsc22d3  2.069058
## 219967     14605   X 140539528 140541107   1580 Tsc22d3  1.362289
## 84835      14229  17  28399095  28401151   2057   Fkbp5  1.596825
## 95516     106878  18  60474191  60475585   1395   Smim3  1.780130
## 89250      57435  17  56103110  56103297    188   Plin4  4.226817
## 174514     53417   7  17030993  17035965   4973   Hif3a  1.681145
## 84836      14229  17  28402609  28402848    240   Fkbp5  1.709063
## 85148      74762  17  29827958  29832399   4442   Mdga1  1.306150
## 126668 100502895   3  88032269  88033152    884 Gm19439  2.455018
## 85004      12575  17  29099344  29100722   1379  Cdkn1a  1.390424
## 29246      19416  11  59963181  59964366   1186   Rasd1  1.732685
## 186521     20887   7 126672870 126673597    728 Sult1a1  1.666580
## 84839      14229  17  28412034  28412124     91   Fkbp5  2.019623
## 102447    226115  19  41063420  41063979    560  Opalin -1.444720
## 19644      12696  10  80170982  80171653    672   Cirbp -1.087081
## 96233      17202  18  66857705  66860472   2768    Mc4r  1.577204
## 84841      14229  17  28428352  28428466    115   Fkbp5  1.998569
##            logCPM        LR        PValue           FDR
## 89251   2.6443083 533.65001 4.538102e-118 7.879961e-113
## 89248   2.2442001 472.55906 8.889999e-105 7.718297e-100
## 89252   1.6642483 320.29747  1.247690e-71  7.221632e-67
## 219969  3.3500243 243.09207  8.327502e-55  3.614968e-50
## 219967  6.1582271 205.72646  1.175647e-46  4.082786e-42
## 84835   5.1146822 189.29529  4.530594e-43  1.311154e-38
## 95516   2.7624054 176.70778  2.536778e-40  6.292658e-36
## 89250  -0.4171852 166.70449  3.881324e-38  8.424413e-34
## 174514  3.2706614 155.59155  1.039736e-35  2.005997e-31
## 84836   1.9322828 144.93730  2.216506e-33  3.848741e-29
## 85148   4.3521796 143.75198  4.025468e-33  6.354385e-29
## 126668  0.8899875 141.03557  1.580388e-32  2.286822e-28
## 85004   4.7712953 133.60940  6.651463e-31  8.884308e-27
## 29246   2.2265890 131.18397  2.256899e-30  2.799199e-26
## 186521  3.1588994 121.08847  3.654654e-28  4.230627e-24
## 84839   1.1119482 107.57110  3.337166e-25  3.621660e-21
## 102447  2.6327847 105.18833  1.110618e-24  1.134398e-20
## 19644   4.7282365 103.29383  2.889517e-24  2.787421e-20
## 96233   2.0361005 101.09557  8.765078e-24  8.010358e-20
## 84841   0.6531217  97.98025  4.225759e-23  3.668804e-19

Alternative Splicing

sp1 <- diffSpliceDGE(fit, contrast=contr, geneid="GeneID", exonid="Start")
## Total number of exons:  173640 
## Total number of genes:  15630 
## Number of genes with 1 exon:  796 
## Mean number of exons in a gene:  11 
## Max number of exons in a gene:  314
topSpliceDGE(sp1, test="Simes", FDR=0.05)
##           GeneID Chr        Symbol NExons      P.Value          FDR
## 78179     545156  16         Kalrn     61 3.668716e-11 5.442173e-07
## 173085 100037396   7 D030047H15Rik      3 1.412498e-10 1.047650e-06
## 15321      20393  10          Sgk1     19 1.602018e-09 7.921447e-06
## 213678     19652   X          Rbm3      8 2.211915e-08 6.269619e-05
## 114243     12064   2          Bdnf      5 2.532590e-08 6.269619e-05
## 150287    231148   5        Ablim2     23 2.535912e-08 6.269619e-05
## 60654      16709  14          Ktn1     42 1.111824e-07 2.356114e-04
## 122598     23856   2         Dido1     17 2.610955e-07 4.841363e-04
## 24849      74309  11         Osbp2     14 3.632241e-07 5.986740e-04
## 126979    319830   3 1500004A13Rik      2 1.342960e-06 1.992147e-03
WriteDataFrameAsTsv(data.frame.to.save=topSpliceDGE(sp1, test="Simes", FDR=1, 
            number=Inf), file.name.path="WTSD5_WTHC5_diff_Splice_edger")

exon plot

plotSpliceDGE(sp1, geneid="Shank3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Kalrn", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Pde4dip", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Pclo", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Fabp7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Hdac7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Per3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Osbp2", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Ablim2", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Ppp2r5c", genecol="Symbol")

plotSpliceDGE(sp1, geneid="D030047H15Rik", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Kalrn", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Sgk1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Dido1", genecol="Symbol")

KOSD5 - KOHC5

Differential expression (QLtest)

contr <- limma::makeContrasts(contrasts="KOSD5 - KOHC5", levels=design)
qlf <- glmLRT(fit, contrast=contr)#glmQLFTest(fit, contrast=contr)
topTags(qlf, n=20, p.value=0.05)
## Coefficient:  -1*KOHC5 1*KOSD5 
##           GeneID Chr     Start       End Length  Symbol     logFC
## 89251      57435  17  56103417  56105501   2085   Plin4  4.154957
## 89248      57435  17  56100591  56102330   1740   Plin4  3.963197
## 89252      57435  17  56106096  56107374   1279   Plin4  3.703562
## 219969     14605   X 140542846 140543185    340 Tsc22d3  2.258351
## 84835      14229  17  28399095  28401151   2057   Fkbp5  1.822751
## 219967     14605   X 140539528 140541107   1580 Tsc22d3  1.458654
## 174514     53417   7  17030993  17035965   4973   Hif3a  1.929814
## 29246      19416  11  59963181  59964366   1186   Rasd1  2.110802
## 89250      57435  17  56103110  56103297    188   Plin4  4.189392
## 95516     106878  18  60474191  60475585   1395   Smim3  1.755609
## 126668 100502895   3  88032269  88033152    884 Gm19439  2.755095
## 84836      14229  17  28402609  28402848    240   Fkbp5  1.804273
## 186521     20887   7 126672870 126673597    728 Sult1a1  1.850060
## 96233      17202  18  66857705  66860472   2768    Mc4r  1.961696
## 85148      74762  17  29827958  29832399   4442   Mdga1  1.278733
## 85004      12575  17  29099344  29100722   1379  Cdkn1a  1.352018
## 84009      71840  17  25471590  25472255    666   Tekt4  3.968303
## 84837      14229  17  28406085  28406270    186   Fkbp5  1.924335
## 9807   100043257   1 150265535 150266066    532 Rbm3-ps -1.149630
## 19644      12696  10  80170982  80171653    672   Cirbp -1.148957
##            logCPM       LR        PValue           FDR
## 89251   2.6443083 546.2779 8.122753e-121 1.410435e-115
## 89248   2.2442001 455.3835 4.858892e-101  4.218490e-96
## 89252   1.6642483 336.8691  3.067190e-75  1.775289e-70
## 219969  3.3500243 279.8074  8.271227e-63  3.590540e-58
## 84835   5.1146822 243.1896  7.929691e-55  2.753823e-50
## 219967  6.1582271 233.7836  8.917448e-53  2.580709e-48
## 174514  3.2706614 200.9858  1.272652e-45  3.156905e-41
## 29246   2.2265890 186.1282  2.225909e-42  4.831335e-38
## 89250  -0.4171852 177.2309  1.950127e-40  3.762446e-36
## 95516   2.7624054 173.3757  1.354926e-39  2.352693e-35
## 126668  0.8899875 171.7133  3.125911e-39  4.934392e-35
## 84836   1.9322828 159.0777  1.799569e-36  2.603976e-32
## 186521  3.1588994 148.6580  3.406455e-34  4.549976e-30
## 96233   2.0361005 143.6157  4.311284e-33  5.347225e-29
## 85148   4.3521796 138.4908  5.691738e-32  6.588756e-28
## 85004   4.7712953 127.4303  1.495646e-29  1.623150e-25
## 84009  -1.0228279 126.1575  2.840277e-29  2.901092e-25
## 84837   1.5768174 123.6899  9.849171e-29  9.501167e-25
## 9807    4.5604543 118.5517  1.312861e-27  1.199817e-23
## 19644   4.7282365 116.4726  3.745089e-27  3.251487e-23

Alternative Splicing

sp1 <- diffSpliceDGE(fit, contrast=contr, geneid="GeneID", exonid="Start")
## Total number of exons:  173640 
## Total number of genes:  15630 
## Number of genes with 1 exon:  796 
## Mean number of exons in a gene:  11 
## Max number of exons in a gene:  314
topSpliceDGE(sp1, test="Simes", FDR=0.05)
##           GeneID Chr        Symbol NExons      P.Value          FDR
## 173085 100037396   7 D030047H15Rik      3 1.427152e-12 2.117038e-08
## 213678     19652   X          Rbm3      8 2.001735e-10 1.484687e-06
## 122598     23856   2         Dido1     17 6.287719e-10 3.109067e-06
## 15321      20393  10          Sgk1     19 1.705489e-09 6.324807e-06
## 78179     545156  16         Kalrn     61 1.347739e-07 3.922097e-04
## 126979    319830   3 1500004A13Rik      2 1.586395e-07 3.922097e-04
## 219971     14605   X       Tsc22d3      5 1.975956e-07 4.187334e-04
## 55935     108143  13          Taf9      7 2.634013e-07 4.884118e-04
## 187121    233908   7           Fus     15 8.310980e-07 1.369834e-03
## 128871     83679   3       Pde4dip     49 1.107124e-06 1.642308e-03
WriteDataFrameAsTsv(data.frame.to.save=topSpliceDGE(sp1, test="Simes", FDR=1, 
            number=Inf), file.name.path="KOSD5_KOHC5_diff_Splice_edger")

exon plot

plotSpliceDGE(sp1, geneid="Shank3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="D030047H15Rik", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Pde4dip", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Pclo", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Fabp7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Hdac7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Per3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Osbp2", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Tsc22d3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Rbm3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Kalrn", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Sgk1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Dido1", genecol="Symbol")

KOHC5 - WTHC5

Differential expression (QLtest)

contr <- limma::makeContrasts(contrasts="KOHC5 - WTHC5", levels=design)
qlf <- glmLRT(fit, contrast=contr)#glmQLFTest(fit, contrast=contr)
topTags(qlf, n=20, p.value=0.05)
## Coefficient:  1*KOHC5 -1*WTHC5 
##       GeneID Chr    Start      End Length  Symbol      logFC      logCPM
## 72408  58234  15 89547626 89549882   2257  Shank3 -9.0760811  5.74935846
## 72409  58234  15 89557735 89560261   2527  Shank3  0.7163855  7.29565025
## 72404  58234  15 89533374 89533454     81  Shank3  1.1667844  2.73493738
## 72654  12805  15 92128128 92128335    208   Cntn1 -3.9086996 -2.18106438
## 72657  12805  15 92220447 92220479     33   Cntn1  1.8140440  0.08011166
## 72414  11434  15 89573832 89574585    754     Acr  1.0215133  1.12615618
## 72402  58234  15 89532095 89532219    125  Shank3  0.7047861  3.09506433
## 72393  58234  15 89503540 89503707    168  Shank3  0.8068976  2.50111953
## 72389  58234  15 89500215 89500418    204  Shank3  0.7701363  2.59030798
## 72421  68708  15 89588188 89588217     30   Rabl2 -1.5991729 -1.27759730
## 72406  58234  15 89543100 89543230    131  Shank3  0.7215936  2.40620387
## 72394  58234  15 89503797 89503913    117  Shank3  0.8909634  1.72778432
## 72400  58234  15 89531343 89531418     76  Shank3  0.9607680  1.51731287
## 72585  66725  15 91814662 91814733     72   Lrrk2 -1.1545411  0.06575861
## 72397  58234  15 89521103 89521373    271  Shank3  0.6465971  3.34820542
## 72403  58234  15 89532329 89532461    133  Shank3  0.6297949  2.28078566
## 72390  58234  15 89501404 89501475     72  Shank3  0.9249444  0.74489213
## 72157  69440  15 89186754 89186819     66 Dennd6b -0.6566132  1.33597099
## 72415  68708  15 89582527 89583430    904   Rabl2 -0.4410104  4.29867874
## 72392  58234  15 89503303 89503454    152  Shank3  0.6778519  1.95573827
##               LR       PValue          FDR
## 72408 2660.68384 0.000000e+00 0.000000e+00
## 72409   71.09526 3.403868e-17 2.955238e-12
## 72404   69.31786 8.380738e-17 4.850771e-12
## 72654   57.58474 3.237270e-14 1.405299e-09
## 72657   45.90491 1.241347e-11 4.310949e-07
## 72414   40.41998 2.048356e-10 5.927943e-06
## 72402   37.01172 1.174216e-09 2.912726e-05
## 72393   33.36129 7.653228e-09 1.484396e-04
## 72389   33.35100 7.693828e-09 1.484396e-04
## 72421   30.24947 3.798950e-08 6.596496e-04
## 72406   29.83175 4.712128e-08 6.871501e-04
## 72394   29.81672 4.748791e-08 6.871501e-04
## 72400   27.68891 1.424775e-07 1.903061e-03
## 72585   24.76712 6.469153e-07 8.023598e-03
## 72397   22.20103 2.455414e-06 2.842387e-02
## 72403   21.84199 2.960504e-06 3.212887e-02
## 72390   21.70289 3.183118e-06 3.251275e-02
## 72157   21.30286 3.921463e-06 3.782904e-02
## 72415   21.12065 4.312550e-06 3.941217e-02
## 72392   20.83469 5.006785e-06 4.346891e-02

Alternative Splicing

sp1 <- diffSpliceDGE(fit, contrast=contr, geneid="GeneID", exonid="Start")
## Total number of exons:  173640 
## Total number of genes:  15630 
## Number of genes with 1 exon:  796 
## Mean number of exons in a gene:  11 
## Max number of exons in a gene:  314
topSpliceDGE(sp1, test="Simes", FDR=0.05)
##       GeneID Chr Symbol NExons      P.Value          FDR
## 72409  58234  15 Shank3     22 0.000000e+00 0.000000e+00
## 70774  27367  15   Rpl3     10 6.015952e-18 4.462032e-14
## 72678  12805  15  Cntn1     26 4.771889e-14 2.359540e-10
WriteDataFrameAsTsv(data.frame.to.save=topSpliceDGE(sp1, test="Simes", FDR=1, 
            number=Inf), file.name.path="KOHC5_WTHC5_diff_Splice_edger")

exon plot

plotSpliceDGE(sp1, geneid="Shank3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Cntn1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Stmn1", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Kif21a", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Rpl3", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Fabp7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Hdac7", genecol="Symbol")

plotSpliceDGE(sp1, geneid="Per3", genecol="Symbol")